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ABSTRACT 

We present numerical studies of the nonlinear, resistive magnetohydrodynamic 
(MHD) evolution of coronal loops. For these simulations we assume that the loops 
carry no net current, as might be expected if the loop had evolved due to vortex flows. 
Furthermore the initial equilibrium is taken to be a cylindrical flux tube with line-tied 
ends. For a given amount of twist in the magnetic field it is well known that once such a 
loop exceeds a critical length it becomes unstable to ideal MHD instabilities. The early 
evolution of these instabilities generates large current concentrations. Firstly we show 
that these current concentrations are consistent with the formation of a current sheet. 
Magnetic reconnection can only occur in the vicinity of these current concentrations and 
we therefore couple the resistivity to the local current density. This has the advantage 
of avoiding resistive diffusion in regions where it should be negligible. We demonstrate 
the importance of this procedure by comparison with simulations based on a uniform 
resistivity. From our numerical experiments we are able to estimate some observational 
signatures for unstable coronal loops. These signatures include: the timescale of the 
loop brightening; the temperature increase; the energy released and the predicted ob- 
servable flow speeds. Finally we discuss to what extent these observational signatures 
are consistent with the properties of transient brightening loops. 
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1. Introduction 

Coronal magnetic loops have been the subject of 
considerable observational and theoretical study. Of 
particular relevanc e here are observations of transient 
brightening loops ( Shimizu et. al. 1994 ). Theoreti- 



cal studies of the stability of such loops have hoped 
to explain some of their observational characteristics. 
In the linear regime the stabilizing effect of line-tied 
boundary conditions has been clearly demonstrated 
( e.g. |Hood fc Priest 1979| ; |Vclli, Einaudi fe Hood 



1990 ). These boundary conditions model the fact that 
on the timescale of ideal MHD instabilities the ends 
of coronal loops can be considered to be frozen into 
a high-density, stationary photosphere. Line-tying al- 
lows coronal loops to have a more twisted magnetic 
field, with consequently more free magnetic energy, 
before the onset of instabilities. However, once a 
critical amount of twist is introduced into the loop 
it does become unstable to ideal MHD modes. The 
non- linear evolution of these instabilities has been the 
subject of a number of recent papers. The ideal evo- 
lution has been studied for a variety of equilibrium 



profiles (Baty & Heyvaerts 1996; Baty 1997; Baty et 



il. 1998). The full resistive evolution has also boon 



studied ( Lioncllo et. al. 1998; Vclli, Lionello & Ein 



ludi 1997 ) but both of these papers assume uniform 
resistivity. In this paper we extend the work of these 
papers by performing detailed non-linear simulations 
from both Eulerian and Lagrangian codes. 

All theoretical studies to date have assumed that 
the initial equilibria are one dimensional, i.e. cylin- 
ders with photospheric line-tying at each end. We 
also make this initial simplification. This has the ad- 
vantage of greatly simplifying the initial equilibrium 
and allows comparison with previous work. We fur- 
ther restrict our attention to loops which carry no 
net current. This class of equilibrium would result 
from slowly applying a coherent vortex flow to the 
flux tube. This has been shown to lead to a coro- 
nal loop of the desired form, with a stable input of 
energy, in numerical experiments (Mikic, Schnack & 



Vaiji Hovcn 1991 ; Van Hovcn, Mok fc Mikic 19135) 



Such an equilibrium, along with the proposed mech- 
anism for its formation, are of course greatly simpli- 
fied, idealisations of the dynamism expected in active 
regions of the solar corona. Here many competing 
effects would be acting simultaneously on each loop. 
Observations show that there is a broad distribution 



are present in active regions. The results here are 
idealised numerical experiments which would be rel- 
evant to the large scale, high energy tail of this dis- 
tribution. This paper is more relevant to micro-flares 
than nano-flares. However, as the distribution of flar- 
ings/brightenings is very broad it should be under- 
stood that for the remainder of this paper when we 
refer to brightening loops we are referring to the high 
energy tail of such events, i.e. compact loop flares or 
micro-flares. 

Several papers have investigated the non-linear 
evolution of a variety of cylindrical equilibria. It is 
important that the relationship between these previ- 
ous studies and our work is made clear. As this paper 
only considers equilibria which carry no net current 
it should be understood that when discussing results 
from these papers we are only referring to those re- 
s ults which are for simila r equilibria . Two papers 
( |Baty fc Heyvaerts 19961 ; |Baty 1997|) have started 
with unstable coronal loops and run codes with large 
implicit steps, steadily increasing viscosity until all 
plasma motion is damped. Eventually this process 
reaches a bifurcated equilibrium state. They have 
found that in this final state the maximum current 
density generated, j m ax, scales with the length of the 
loop. In this paper we resolve the full time-dependent 
solution and find no evidence that j ma x scales with 
the loop length. We find that for simulations without 
resistivity j m ax is always the largest current possible 
for the given resolution. This is determined by the 
numerical mesh size in the region of jmax, not the 
loop length, and is true for both the Lagrangian and 
Eulerian simulations. This does not contradict the re- 
sults about equilibrium current densities as the j max 
we observe is not part of an equilibrium. Indeed, we 
observe these large current concentrations at a time 
when there are velocities (resulting from the insta- 
bility) of the order of the local Alfven speed. Such 
features would have been damped by the numerical 
procedure adopted in searching for equilibrium solu- 
tions. The code used to find these equilibrium solu- 



tions has also been run with a fixed viscosity ( Baty et 



of energy release events ( Shimizu et. al. 1994 ) which 



al. 1998 ). This paper found that the ideal instability, 
and consequently jmax-, saturated before the current 
density reached the grid scale length. This is in con- 
tradiction with the results presented here in which 
there is no evidence for ideal MHD saturation of the 
instability. Furthermore, the resistive phase which is 
triggered by the large current densities in our simu- 
lations also prevents the loop from ever reaching the 
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ideal MHD bifurcated equilibria described earlier. In 
this regard we are in agreement with other non-linear 
simulations (Lioncllo et. al. 199S; Vclli, Lionello & 



Eiraudi 1997) which show the same collapse to grid 
scale lengths resulting from the ideal MHD instabil- 
ity, i.e. no signs of saturations of the ideal mode. We 
point out here that while it is often assumed that cur- 
rent sheets form as a result of ideal MHD instabilities 
this has never actually been proven. Some papers 
(Baty et. al 1998| ) actually contradict this belief 



while the only tw o papers which seem to lend weight 
to this argumen t ( Lioncllo et . al. 1998|; |Velli, Lioncllo 



& Einaudi 1997 ) have only published results for a sin- 
gle resolution. At this resolution the current density 
in the current concentration is only a few times that 
present in the initial equilibrium. Bearing in mind the 
lack of consensus on this issue we deal with the for- 
mation of current sheets in some detail in this paper 
before treating the resistive evolution. The use of a 
Lagrangian code in this paper allows the formation of 
current densities which are two orders of magnitude 
larger than has been possible in previous studies. We 
also present detailed tests of the scaling of the current 
density with grid size for both Lagrangian and Eule- 
rian simulation. These tests combine to give the most 
convincing evidence to date that current sheets do in- 
deed form as a result of ideal MHD instabilities in the 
corona. In this paper what is meant by a current sheet 
is that unless some dissipative processes is introduced 
into the system the current density will continue to 
collapse to shorter scale lengths without limit. In re- 
ality it is resistive effects which stop this collapse. It 
is the treatment of resistivity, the greater number of 
scalings on different grids and tests with a Lagrangian 
code which distinguishes our work from previous pub- 
lications. It is also these differences which allows us to 
make the first quantitative predictions about unstable 
coronal loops. 

While the papers discussed in the previous para- 
graph have treated a variety of equilibria here we 
study only equilibria which carry no net current. We 
concentrate on a single representative equilibrium and 
perform a variety of detailed numerical experiments. 
In these we abandon the common practice of assum- 
ing a uniform resistivity and instead couple the re- 
sistivity to the current. In this model the resistivity 
is only present in regions where the current density 
exceeds some critical value. In this way resistivity is 
only applied in those regions where reconnection is 
allowed. This procedure allows us to get closer to a 



converged answer for real coronal values and gives us 
increased confidence in the quantitative predictions of 
observational signatures which we make from our sim- 
ulations. The evolution is largely independent of the 
form of resistivity used provided that it is localised 
in the current sheet. This point has been postulated 
in other papers flVclli, Lioncllo fc Einaudi 1997 ) but 
is confirmed for the first time in our simulations. We 
also include a comparison with simulations in which 
the resistivity is assumed constant and highlight the 
differences between the two cases. 

The paper is organised as follows. In §2 we de- 
scribe the physical model and explain our motivation 
for the choice of resistivity. §3 contains the details of 
the equilibrium we have chosen along with its linear 
stability properties. The formation of current sheets 
has been investigated numerically and the results of 
this study are included in §4. The full non-linear re- 
sistive evolution of the loop is outlined in §5. Finally 
the conclusions which we can draw from our stud- 
ies, including predictions of observational signatures 
of unstable loops, are presented in §6. 

2. Physical Model 

In this paper we represent coronal loops as initially 
cylindrical tubes. These tubes have each end tied into 
the photosphere. Photosphcric line-tying is modelled 
by imposing zero velocity at the ends of the loop along 
with suitable symmetry properties on the components 
of the magnetic field. We ensure that the tube, of 
length L z , is sufficiently long that it is unstable to 
ideal MHD modes. The solution domain is Cartesian 
with transverse size (L x , L y ). The instability remains 
localised inside a specific radius so that provided (L x , 
L y ) are large enough the choice of boundary condi- 
tions in (x, y) makes no difference to the evolution of 
the loop. 

The evolution of the coronal loop is modelled by 
the resistive MHD equations. In dimensionlcss form 
these are 
dp 
dt 

dt 
dB 

~dt 

de 
dt 
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locity , P is the thermal pressure, e = P/(F— 1) is the 
internal energy density (T = 5/3 is the specific heat 
ratio), p is the mass density and r\ is the resistivity. 
For all simulations we assume a plasma (3 of 0.01. This 
model ignores thermal conduction, radiation, grav- 
ity and heating terms other than Ohmic. Ignoring 
the transport terms is justified because of the short 
timescales involved in these simulations. The validity 
of the neglect of gravity is of course dependent on the 
actual length of the loop being considered. With pres- 
sure scale heights of ~ 100 Mm in the corona this is 
a good approximation for transient brightening loops 
but clearly less valid for large quiescent loops. 

The choice of the functional form for the resistivity 
r\ in the simulations presented in this paper is partic- 
ularly important. Normally r\ is treated as a constant 
chosen on numerical grounds such that the explicitly 
included resistive diffusion exceeds numerical diffu- 
sion. In this paper rj is chosen so that it is only in- 
cluded in those regions where it is needed. The clas- 
sical resistivity of coronal plasmas is negligible on the 
timescale of MHD instabilities, i.e. the timescale of 
interest here, everywhere except in regions of intense 
current concentration. Hence we choose a resistivity 
given by 

r, = r) O MAX(0,\j\-j crit ) (5) 

In this formula rj = if \j\ < j C rit where j cr u is the 
critical value of current density needed before the re- 
sistivity is turned on. Thus while 770 is still chosen 
on numerical grounds the resistive effects are only 
applied to those regions in which they are actually 
needed. The point here is essentially that if we were 
to take a uniform resistivity we would be diffusing 
the magnetic field in regions were resistivity should 
be negligible. It is true that locally the effect of uni- 
form resistivity would be largest in the current sheet 
and in that region the difference due to the choice 
of uniform resistivity over equation || would be small. 
However, the small resistive effects which are then ap- 
plied over a much larger total volume do change the 
nature of the final solution. This will be discussed 
in §5 where a comparison of the two approaches is 
presented. 

If intense current concentrations form as a result of 
the ideal MHD instability then the electron fluid flow 
speed U e will become large. These current concentra- 
tions lead to localised heating of the electrons and on 
the timescales of interest here these are not in thermal 
equilibrium with the ions. Thus the electron temper- 
ature will exceed the ion's and once U e exceeds the 



local ion sound speed the current concentration will 
drive ion-acoustic turbulence. The fluctuating elec- 
tric field of this turbulence causes electron scattering 
which is manifested at the fluid level as an increase in 
the local resistivity ( Rosner et. al. 1978| ). By rapidly 
creating a current sheet the ideal MHD instability 
therefore creates all of the conditions necessary for 
the onset of ion-acoustic turbulence. These are that 
the electron temperature exceeds the ion temperature 
and that the electron flow speed exceeds the local ion 
sound speed. While these conditions are satisfied for 
the current sheets driven by ideal MHD instabilities 
it must be noted that in other circumstances in which 
current sheets form this may not be the case. If the 
electron flow speed continues to be driven, as is the 
case here, this enhanced resistivity increases. Equa- 
tion H would then be a suitable macroscopic param- 
eterisation of the sub-grid scale turbulence. For a 
coronal loop with number density 10 16 to -3 ; ion ther- 



mal speed 1.3 x 10 5 



width 10 m and magnetic 



field strength 100 G the normalised critical current 
density required for the onset of such turbulence is 
jcrit ~ 3000. These are order of magnitude estimates 
of coronal values appropriate for transient brighten- 
ing loops ( Shimizu 1996| ). Unless otherwise stated 
these characteristic values for density etc. are the 
ones used throughout the paper when estimating real, 
i.e. unnormalised, values. Note however that the tem- 
perature assumed in calculating the thermal speed is 
2 x 10 6 K. This is taken as an average temperature 
of an active region and is meant as an estimate of 
the loop temperature before the brightening occurs. 
Of course the codes use normalised variables so scal- 
ing to other coronal values is a trivial matter and 
does not affect the qualitative features of these sim- 
ulations. This estimate of j cr u should be compared 
with the equilibrium current density which has a max- 
imum value of 4.5 in these units. When ion acoustic 
turbulence is active the characteristic electron scat- 
tering time is approximately the ion plasma period. 
This gives an increase in the plasma resistivity by 
a factor of up to 10 6 . If one uses the Spitzer for- 
mula for classical resistivity one concludes that the 
normalised resistivity is ~ 10~ 12 . However in cur- 
rent sheets which trigger ion-acoustic turbulence 770 
in equation (5) is ~ 10~ 6 . It should also be pointed 
out that ion-acoustic turbulence is not the only pos- 
sible source of enhanced resistivity. For example the 
initial rapid heating in the current sheets of unsta- 
ble loops has been shown to trigger Langmuir wave 
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enhanced resistivity (Takakura 1991). While there 
is theoretical evidence that the mechanism described 
above would lead to enhanced resistivity there is no 
direct evidence that such levels of turbulence are in- 
deed generated and sustained on the required scale 
in current sheets in the solar corona. As a result the 
arguments in support of the use of Equation || as the 
correct physical form of resistivity, as opposed to be- 
ing simply numerically appropriate as discussed in the 
previous paragraph, must remain speculative at this 
stage. 

3. Equilibrium and Linear Stability 

Throughout this paper we limit our attention to 
just one equilibrium. This has allowed us to perform 
a detailed set of numerical experiments, over a range 
of grid sizes, using different codes. The equilibrium 
is taken to be a force-free cylinder which carries no 
net current. The precise form is given in terms of the 
axial current density in the loop, j z , by 



(6) 



The normalisation used is that the loop is confined 
within a radius r = 1 so that the free parameters a 
and b are found from the conditions = and 

Jo Tj z {r)dr = 0. These guarantee that the equilib- 
rium does not contain a surface current at r = 1 and 
that the total current in the loop is zero. The actual 
values are a — 5/(6^) and b = 1/V&- The poloidal 
component of the magnetic field, Bg, is then found 
from Ampere's law and B z from the force-free condi- 
tion, 

Bl-Bl-2f (?f)dr' (7) 



Bi 



This equilibrium definition has two free parameters, 
Bq and jo- These are chosen to be 1.0 and 4.3 re- 
spectively. With this choice the equilibrium profiles 
are shown in Figure Al. Note that outside r = 1 the 



magnetic field is purely axial. This choice of equilib- 
rium is similar to equilibria with no net current stud- 



ied elsewhere (Velli, Lionello & Einaudi 1997; Baty 



t. al. 1998; lionello et. al. 1998). These works 



all showed that once the critical length for the on- 
set of instability has been exceeded the growth rate 
increases rapidly, eventually reaching a growth rate 
which is almost independent of length. Our choice 
of L z = 10 for most of our simulations puts the 



equilibrium in this region. We therefore do not self- 
consistently follow the loop as, for example, its length 
increases from a stable length up to the length con- 
sidered here. We simply assume that it is initialised 
at a length which is already unstable at the saturated 
growth rate. This is necessary numerically so that 
the instability grows sufficiently quickly to be simu- 
lated in a reasonable time. We have performed some 
simulations for shorter initial lengths and found that 
our predicted observational signatures are insensitive 
to this choice. As far as the ideal phase of our simu- 
lations are concerned we show broad agreement with 
the results in earlier works. This shows that these 
results are not critically sensitive to the details of the 
equilibrium and is a useful confirmation of the results 
from a completely different set of codes. This paper 
differs from those earlier publications in performing 
more detailed scaling tests on the formation of the 
current sheet and more detailed simulations of the re- 
sistive phase. The choice of equation || for the resis- 
tivity is a particularly important difference between 
this work and previous publications. 

In the first instance we studied the linear ideal 
MHD stability of our model equilibrium. The MHD 
equations (0)-(|J) for zero resistivity are linearised, 
and a time dependence oc exp(7t) assumed to ob- 
tain the equations for the linear normal mode spec- 
trum in ideal MHD. These normal modes are then 
calculated using a bicubic finite clement code similar 



to that described in Van dcr Linden & Hood 1998 



The qualitative features of the stability of this equi- 
librium are the same as those for previously studied 
equilibrium profiles, e.g. Lionello ct. al. 1998. Un- 



less otherwise stated the length of the loop in the re- 
mainder of this paper will be fixed at L z — 10. This 
point is chosen some distance from the marginal sta- 
bility point, so that the growth rate of the fundamen- 
tal mode (7 » 0.302) is nearly equal to its maximal 
value. This allows the instability to develop suffi- 
ciently rapidly to save computational overhead. At 
this loop length, two more instabilities exist: the first 
overtone has a growth rate of 0.197, while the sec- 
ond has growth rate 0.037. The latter is certainly 
not expected to have any significant contribution to 
the development of the instability due to its small 
growth rate, while the former might contribute, but 
in practice does not show up in the time-dependent 
simulations presented in this paper (which is mainly 
due to the choice of initial perturbation) . 

To make a clear distinction between the growth of 
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the linear instability and the effects of non-linearity 
on its evolution, it is very useful to compare both 
linear and non-linear time-evolution of the same ini- 
tial perturbation. A good approximation of the linear 
time-evolution may be obtained from the spectrum of 
normal modes by writing the initial perturbation as 
a linear combination of the full spectrum of normal 
modes. A quantitative comparison of the linear ver- 
sus non-linear evolution is made at t = 5 in Figure 
A.2, where the perturbed current density is plotted 
along a line through the centre of the loop parallel 
to the photospheric line-tied ends from the linear and 
both non-linear simulations. It is clear that at this 
time, the evolutio n of the instability is still linear in 
nature. In Figure A3 a similar comparison is made 
at the later time of t = 10. At the early time, i.e. 
t = 5, the linear and non- linear results differ in two 
noticeable ways. Firstly the peaks from the non-linear 
analysis are shifted to the right due to the central 
plasma column moving in the x direction. Figure A4 
shows v x along the same a;-axis at t = 5 from the 
linear code. From this one can see that the regions 
with the steepest gradients are located in the region 
0.6 < |x| < 0.8. These are regions of compression for 
x > and expansion for x < 0. They account for 
the asymmetry seen in both non-linear results. By 
t = 10 the general structure of the linear mode is 
still evident but non-linear effects are becoming sig- 
nificant. Most prominent is the formation of a large 
current density at x = 0.8. Note that Figure A3 is 
truncated at \ji\ — 10 and that for the Lagrangian 



code |ji 



= 31 at t = 10 



4. The Formation of Current Sheets 

The fully non-linear evolution of the instability is 
followed using two time-dependent nonlinear MHD 
codes. The first is an ideal MHD Lagrangian code 
that follows the initial phase and the formation of the 
current sheet, the second a resistive MHD Eulerian 
code that allows the resulting reconnection to be fol- 
lowj od and the later phases of the evolution examined. 
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code of Longbottom et. al. 1998 and is described 



in the appendix. It has two main features relevant 
to this study. The code solves the ideal MHD equa- 
tions (there is no dissipation due to resistivity or vis- 
cosity). The grid on which the equations are solved 
moves with the fluid and thus in this case, where the 
inner region of the loop is forced out against the near 



stationary outer potential field, more grid points will 
accumulate at the regions where large gradients form. 
These two features together ensure that the current 
structures resulting from the instability remain highly 
resolved. 

The loop is confined inside r = 1 and centred in 
the (x, y) computational plane which is L x — 6 by 
L y — 6. The simulations were carried out on grids of 
61 x 61 x 21, 91 x 91 x 31 and 151 x 151 x 101 points 
in the (x, y, z) directions. The [x, y) grid is uniform 
inside —1.1 < x, y < 1.1 and only stretched outside 
this central core, with the initial grid spacing (at t = 
0) within -1.1 < x,y < 1.1 being 0.05, 0.033 and 
0.02 respectively. The z grid is uniform with —5 < 
z < 5. Simulations are started with a small velocity 
perturbation (v max = 0.01) whose structure is taken 
to approximate that found from the linear analysis. 
The results below are shown for 151 x 151 x 101 grid 
points. 

As described in the previous section the initial evo- 
lution agrees well with that predicted by linear theory, 
the linear eigenfunction and growth rate being repro- 
duced. However, from t = 5 onwards a helical cur- 
rent structure grows as a result of the inner twisted 
magnetic field being forced against the surrounding 
potential field by the developing kink mode. This be- 
haviour can be seen in Figure AE which shows the 
perturbed current ( \j\ \ ) plotted along the x-axis at 
y = z — as a function of time. Here the maximum 
value of current plotted has been truncated at j\ = 50 
so that the details at earlier times can be seen. The 
actual maximum value of current in the current sheet 
at t = 11.5 is ji = 766. Even at late times the rem- 
nant of the linear mode are still visible in the central 
part of the loop. The formation of the current struc- 
ture at the rational surface can be clearly seen. The 
maximum current after t = 10 scales faster than n x , 
where n x is the number of gridpoints in the x direc- 
tion in the Lagrangian code. No sign of saturation 
of the current is seen. This behaviour is indicative of 



the f ormation of current sheets ( Longbottom et. al 



1998| ) 



A surface plot of the total current in the (x, y) 



plane at the loop apex is shown in Figure A£ . Again 
the current has been truncated at j = 50 so that both 
the shape of the current concentration and the inner 
structure is visible. The global structure is that of 
a helix wrapped around the central loop column and 
is essentially identical to that found by the Eulerian 
code, Figure A7 a. The half- width of the current sheet 
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at t — 11.5 is approximately 0.0025. This would mean 
that for the Eulerian code to resolve this current at 
this time there would need to be ~ 3500 grid points 
in both the x and the y directions. 

These results have also been confirmed with the 
Eulerian code described in detail in the next section. 
The important point as far as the current sheet for- 
mation is concerned is that running in the ideal MHD 
model, i.e. with r\ = 0, the current generated in the 
current sheet scaled as 1/dx, where dx is the grid 
resolution in the current sheet, and showed no signs 
of saturating as the resolution was increased. Note 
that the scaling of the maximum current is different 
than in the Lagrangian code as the Lagrangian code 
moves points into the region where the current sheet 
is formed. For the Eulerian code the highest current 
obtained in the current sheet was ~ 30. This be- 
haviour is consistent with previous work in this area 
( iLioncllo ct. al. 1998| ; |Vclli, Lioncllo fe Einaudi 1997p 
although by performing scalings on different grid res- 
olutions we are able to verify the correct scaling with 
mesh size which is required if these are indeed cur- 
rent sheets. Furthermore the current densities gener- 
ated in the Lagrangian code are two orders of magni- 
tude larger than has been possible with the Eulerian 
codes used in previous studies. The results in this 
paper contradict some simulations which found non- 
linear saturation of the instability with a purely ideal 
MHD description flBaty et. al. 1998|). This discrep- 



ancy between different codes has been noted before 
(|Baty et. al. 1998|; [Lioncllo ct. al. 1998j) and has 



been attributed to the different treatment of small 
scales, with associated numerical dissipation, in the 
codes used. Here we have used two non-linear codes, 
each with distinctively different properties from ear- 
lier codes, and found no saturation. This is the first 
time a Lagrangian code has been for such simulations 
and that tests on the scaling of the current density 
with dx have been presented. These provide strong 
evidence that these instabilities do not saturate while 
still described by ideal MHD and that the gener- 
ated current concentrations are current sheets. We 
have also performed simulations of loops with lengths 
L z = 5.5 finding the same magnitude current sheet 
evolving in both cases. Previous studies of bifurcated 
equilibria ( Baty fc Heyvaerts 1996 ; Baty 1997 ) have 
found that j ma x scales linearly with L z . Our sim- 
ulations show that this is not true of the instability 
driven current sheet which forms while the loop is still 
in a highly dynamic, non-equilibrium state. 



5. The Resistive Phase 

The Lagrangian code is only valid for ideal MHD. 
To follow the resistive evolution of the loop an Eule- 
rian code is used. The ideal MHD part of this c ode is 
based on the MH3D code ([Lucek fc Bell 1996|) writ- 
ten at Imperial College, London. The code uses a 
stretched Eulerian grid. Variables on this fixed grid 
are updated using a split Lagrangian, Eulerian remap 
technique, with the advection handled by a second or- 
der Van Leer upwind scheme ( Youngs 1982 ). The 
code maintains V.B = by using the constrained 
transport model for magnetic flux advection (Evans & 



Hawley 1988). During the Lagrangian phase of each 



timestep an artificial viscosity term is added to equa- 
tion |^. This is applied as a viscous pressure at cell 



boundaries for cells which are being compressed (van 



Neumann & Richtmyer 1950). This viscosity results 
purely from compressive effects, i.e. shear viscosity 
is not included. For comparison with papers which 
add a viscous term of the form ^/?V 2 v to equation]^ 
the formula we use is approximately equivalent to this 
form with v = 10~ 3 . We have confirmed this with di- 
rect comparison of the two forms on low resolution 
runs. It should also be pointed out that in our units 
v = 10 -3 is the correct order of magnitude for a tran- 
sient brightening loop. The viscous heating in these 
simulations is significant in the overall energy balance 
and so must be included in the energy equation. 

The results presented here are from runs with a 
161 3 Cartesian grid. This is the resolution used in 
the largest set of numerical experiments and there- 
fore constitutes our largest consistent data set. Some 
higher resolution tests have been performed to test 
convergence. These had a 221 x 221 x 101 grid in 
(x, y, z) with the grid stretched to give twice the reso- 
lution in (x,y) of the 161 3 experiments. Results from 
this resolution will be called the high resolution re- 
sults in the remainder of this paper. Unless explic- 
itly stated it should be assumed that results are from 
the 161 3 runs. The loop is confined inside r = 1 
and centred in the (x, y) computational plane which 
is L x — 6 by L y — 6. The (x,y) grid is uniform 
inside —1.1 < x,y < 1.1 and only stretched outside 
this central core. The ratio of minimum to maxi- 
mum grid spacing was 3.2. The z grid was uniform 
with the coordinate range —5 < z < 5. Simula- 
tions are started with a small velocity perturbation 
{v-max = 0.01) whose structure is taken to approxi- 
mate that found from the linear analysis. 
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Figure A7 shows iso-surfaces of the magnitude of 
the current density at two different times. The sur- 
faces are iso-surfaces of |j| = 3. As the maximum 
current in the initial equilibrium is 4.5 these sur- 
faces show the perturbed central column as well as 
the current sheet. While current density is not ex- 
perimentally observable these figures are presented 
as they give the clearest picture of the physical pro- 
cesses present in this unstable loop. These results 
are from a simulation with 770 = 10~ 3 and j C rit = 5. 
The value of 770 is the smallest value we can use in 
this code and still guarantee that the deliberately 
included, numerically controlled resistivity is larger 
than the numerical resistivity inherent in the differ- 
ence scheme. j cr it is fixed to the largest value that 
allows us to fully resolve the current sheet. The non- 
linear feedback through resistivity of the form given 
in equation then restricts the maximum current in 
our current sheets to about 10. Figure A7(a) shows 
the iso-surface at t = 10 with the central column per- 
turbed into the characteristic helical m — 1 mode. 



Wr apped around this central column is the current 



sheet which is formed at the place where the pitch of 
the instability matches the pitch of the magnetic field 
(see Baty & Heyvaerts 1996 for a detailed discussion 
of this process) . At this time the current in this outer 
current sheet has just reached 5. Therefore up to this 
point equation || has been setting rj = everywhere 
and the code has been solving the ideal MHD equa- 
tions. Beyond t = 10 the current in the current sheet 
continues to increase and equation |5| 'turns on' the 
resistivity and the code automatically begins solving 
the resistive MHD equations but with the resistiv- 
ity only present in the outer current sheet. Figure 
A7(b) shows the current iso-surface at t = 15. More 
of the central column has been moved out, by the 
ideal MHD instability, into the region of the current 
sheet. In this region reconnection is allowed so that 
the twist in field lines can be removed, or equivalently 
the current dissipated. This process continues until 
at t — 20 sufficient current has been dissipated that 
no region has \j\ > 5 and the resistivity 'turns off. 
In summary Figure A7 shows that the ideal MHD in- 
stability drives the current in the loop out into the 
current sheet were it is dissipated. 

The above experiment has been run with two dif- 
ferent plasma density profiles. In both cases the pres- 
sure is uniform, as required by the force- free con- 
dition. In one set of tests the plasma density was 
taken to be uniform across the whole computational 



domain. In the other the density profile was taken to 
be p — 0.45(l+cos(7rr)) + 0.1 for r < 1 and p = 0.1 for 
r > 1. This second choices makes the average density 
inside the loop 3.67 times that of the surrounding 
coronal plasma and was motivated by observations 
that the density inside a brightening loop exceeds 
that of the surrounding coronal plasma. These ex- 
periments showed that the evolution of the unstable 
loop is insensitive to the choice of density profile. It 
should be noted that while the second density profile 
docs imply a drop in temperature inside the loop this 
is unimportant for MHD simulations. It is the pres- 
sure which exerts a force and the temperature does 
not appear in equations (l)-(4). Changing the den- 
sity merely changes the timescales involved. What is 
clear from the second density profile is that in the fi- 
nal state, i.e. at t = 20, the density enhanced region 
still lies within the same bounding radius. In other 
words the loop is not destroyed by the instability but 
the twisted magnetic field lines are straightened out 
This is in agreement with previous studies (Lioncllo 



et. al. 1998) 



Figure A8 shows two different iso-surfaces of en- 
ergy density taken at the same time, t = 20. Fig- 
ure [A^(a) shows regions which have been heated to 3 
times the initial background value. This shows that 
the loop has brightened along its whole length. Fig- 
ure A8(b) shows the higher energy components which 
have been heated to 6 times the background value. 
In this Figure the central region has been heated by 
Ohmic dissipation in the current sheet while the ends 
have been heated by the viscosity in the code. The 
rings at each end of Figure |A|(b) are where the loop is 
tied into the photosphere. The very dynamic nature 
of the loop at t = 20, see below, causes viscous stresses 
at the ends where the velocity is forced to be zero by 
the photospheric boundary conditions. Once the cur- 
rent density drops below j C rit everywhere the resistive 
phase is over and the code reverts to ideal MHD (al- 
though viscosity is still present). No simulations have 
been performed beyond this time as from that time 
onwards the timescale of interest is the timescale for 
thermal conduction. This is too long to be studied in 
3D with this kind of resolution. 

At t = 20 there are very large flows set up due to 
the ideal instability and field line reconnection. The 
peak value is ~ 0.8 Va, where Va is the Alfven speed. 
However, observations of flows in the corona are based 
on Doppler shift measurements. Such measurements 
include averaging over exposure times, pixel sizes and 
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line of sight effects. To estimate the importance of 
these effects we have taken a simple density weighted 
average of one of the transverse components of ve- 
locity, i.e. v x or v y , over an area of approximately 
1.5 Mm x 1.5 Mm, an exposure time of 10 seconds 
and along a line of sight. This averaging time and 
area are typical of solar observations. The flow struc- 



ture after this averaging is shown in Figure A9 . In this 



figure the photospheric footpoints are at ±5Mm and 
the points are the centres of our pixels. This sort of 
simple averaging can only be taken as an estimate of 
the kind of velocities which could be observable exper- 
imentally Issues such as ionisation population levels, 
temperature dependence of the weighting fun ctio n are 
beyond the scope of the current work. Figure A£ does 
however show that the very large flow speeds present 
in our simulations would not be directly observable. 
These simulations therefore suggest that after a loop 
with parameters typical of a brightening loop has gone 
unstable, flows of ~ 40km s^ 1 should be observable. 
The averaging implicit in these measurements would 
however mask the real loop plasma flows which are 
highly localised and as large as 1500km s^ 1 . 

The energy released from this instability is 54% of 
the available magnetic energy. This is split almost 
equally between Ohmic heating, kinetic energy and 
viscous heating. The available energy is defined as the 
energy stored in the Bg component of the equilibrium 
magnetic field. For this equilibrium the free magnetic 
energy for the coronal values in §2 is 9.5 x 10 28 ergs. 
For typical brightening loop values the resistive phase 
lasts about 5 seconds. 

For comparison we have repeated the above simu- 
lation with a uniform resistivity. This is turned on at 
t = 10 and the simulation is stopped at t = 20. In this 
way the resistivity is applied for the same length of 
time as above. For this run we took 77 = 10~ 3 so that 
this too is consistent with the value used above. We 
find that for this uniform resistivity model 62% of the 
available energy is released. This is sufficiently close 
to the value found from using equation ^| that the dif- 
ference can be ignored. However, the kinetic energy 
generated with a uniform resistivity is approximately 
half that of the value from using equation ^| and the 
total Ohmic heating is three times larger. The peak 
flows are less than half those shown in Figure pvTj| . 
At present all simulations are forced to use resistivity 
which is unphysically large. These tests show that 
applying such a large resistivity uniformly over the 
computational domain, instead of localising it to just 



those regions where it should have an effect, over- 
estimates the Ohmic heating and under-estimates the 
kinetic energy in the final dynamic state. Neither of 
these points is surprising as including resistivity ev- 
erywhere will clearly increase the overall Ohmic dis- 
sipation and smooth the fields driving the instability. 
What is important here is that these effects have now 
been quantified for the values of resistivity typically 
used in large scale numerical simulations. Using a 
uniform resistivity does release the same amount of 
magnetic energy as using equation || but splits it be- 
tween Ohmic heating and kinetic energy in a very 
different way. It is worth noting that runs with zero 
resistivity also release the same amount of energy but 
most of this energy is simply lost from the system. 
In this case when the current density scale length 
reaches the grid spacing numerical diffusion dissipates 
the current. Thus getting the correct amount of to- 
tal magnetic energy released in a simulation is not of 
itself a useful indicator that the resistive effects have 
been correctly modelled. 

All of the results above are from 161 3 stretched 
grids with 770 = 10~ 3 and j C rit = 5. By running the 
same simulations on 81 3 and 12 1 3 grids we have con- 
firmed that, for these values of 770 and jcrit, these 
results are the correct, converged solutions. We have 
also conducted three higher resolution simulations on 
a stretched 221 x 221 x 101 grid. The first was a 
purely ideal MHD simulation. This confirmed that 
the current density in the current sheet scales as 1 / dx 
demonstrating that we have no evidence of the current 
density saturating consistent with the results from 
the Lagrangian code. The second simulation repeated 
previous resistive runs with rjo = 10~ 3 and j cr u = 5. 
This confirmed the accuracy and convergence of these 
results. The last high resolution run repeated this 
last simulation but with j C r%t = 10. It is only at this 
higher resolution that enough grid points are present 
in the current sheet for energy conservation to be ac- 
ceptable, i.e. the Ohmic heating is much larger than 
the energy loss through numerical diffusion, for this 
value of j cr it. The interesting point here is that the 
results with j cr u = 10 are not significantly different, 
i.e. none of the observational signatures change, from 
those with j C rit = 5. Similarly, increasing 770 to 10 -2 
does not alter our observational predictions. Hence 
over the range of dimensionless parameters resolvable 
by the simulations we have performed we find that 
the predicted observational signatures are insensitive 
to the choice of 770 and jcrit- However, the use of 
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equation [5| for coronal values assuming that turbu- 
lence enhanced resistivity (as discussed in §2) was ac- 
tive would require 100 times the resolution we have 
used. If such enhanced resistivity were absent then of 
course much higher resolution still would be required. 
Such grid sizes are of course impossible at present. 
Our predicted signatures are therefore based on the 
assumption that this independence on rjQ and j cr n re- 
mains true up to coronal values. 

The results above are for a loop with L z = 10. 
We have also checked this result for L z = 5.5 which 
has approximately half the growth rate and is closer 
to the marginal stability length. The resistive evo- 
lution of this loop showed the same behaviour as for 
L z = 10. The maximum current in the current sheet 
was the same; the timescale for the resistive phase 
was 10 Alfven transit times as before and the en- 
ergy released was the same fraction of the available 
free energy. The only difference was that as a re- 
sult of the lower growth rate it took longer to reach 
the stage where the current density in the current 
sheet triggered the resistivity. A final set of tests has 
also been performed in which the resistivity was set 
to a constant value of 10 -3 in any computational cell 
with current density larger that j crit . These also gave 
the same set of observational signatures as presented 
above verifying that the results are also insensitive to 
the functional form chosen for the resistivity provided 
that it is only present in the current sheet. 

6. Discussion 

Our aim in this work has been to perform non- 
linear numerical solutions of unstable coronal loops 
using resistive MHD. In these simulations we have 
concentrated on a single equilibrium. We have con- 
firmed some of the results of previous papers using 
different initial conditions and different numerical ap- 
proaches. Where there has been a lack of consensus 
in previous studies we have clarified these issues by 
supplying detailed numerical results. These include 
the formation of current sheets; the basic structure 
of those current sheets and the general features of 
the resistive evolution. There are however important 
quantitative differences between this paper and pre- 
vious work. Most of these differences stem from only 
applying resistivity locally in current sheets. Before 
progressing to a description of the observational sig- 
natures it is worth highlighting what those differences 
are. 



1. We have run the non- linear codes on a range 
of grids. These confirm that the current gener- 
ated near the resonance surface scales as 1/dx 
for the Eulerian code and faster than 1/dx 2 
for the Lagrangian code. The current densities 
found in the Lagrangian code are two orders of 
magnitude large than those found in previous 
studies. This combined with the above scalings 
with grid size provide convincing evidence that 
the current density would continue to collapse 
to smaller and smaller scales unless stopped by 
some dissipative process, i.e. they are current 
sheets. 

2. Coronal simulations require the use of a resis- 
tivity which is larger than the real coronal value 
in order to limit the current densities formed in 
such simulations. We adopt this procedure but 
only apply resistivity where it is needed by us- 
ing equation |^ to localise resistive effects to re- 
gions with high current densities. If the current 
density exceeds 3000 (in our normalised units) 
then equation || may be viewed as a parameteri- 
sation of the effects of sub-grid scale turbulence 
enhanced resistivity based on theories of ion- 
acoustic turbulence. While there is theoretical 
evidence that this may be true a direct proof of 
this is beyond the scope of this paper. 

3. We have quantified the discrepancies between 
using equation ^| and taking rj to be uniform 
over the computational domain. 

4. For all of the tests we have run, the observable 
properties of the loop are independent of the 
choice of j cr u and 770 in equation [|. As far as 
we can determine it is the rate at which mag- 
netic flux is moved into the reconnection region 
which is important. This is determined by the 
correct resolution of the ideal MHD instability. 
We have also shown that the results are insen- 
sitive to the function form of equation 5. 

The simulations we have performed are only for one 
equilibrium. This equilibrium carries no net current. 
Such an equilibrium might evolve due to long scale 
length correlated twisting of the flux tube either in 
its rise through the convection zone or through pho- 
tospheric vortex flows once it has emerged into the 
corona. We further assume that the equilibrium is 
force-free and that at the start of our simulation the 
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loop is unstable. The loop could have become unsta- 
ble from increased twisting of the magnetic field due 
to photospheric motion or from rising higher into the 
corona and hence increasing in length. Provided these 
conditions are satisfied, or at least a reasonable ap- 
proximation, then the observational signature of the 
ensuing instability are those set out below. Where 
results are presented in unnormalised units it is as- 
sumed that they relate to the typical brightening loop 
values listed in §2. While we have shown that our re- 
sults are independent of tjq and j cr u over a large range 
of values it should be noted that these observational 
signatures are only valid if we assume that this re- 
mains true when extrapolating to real coronal values. 
Implicit in this is of course the assumption that the 
microscopic details of magnetic field diffusion and re- 
connection can be adequately parameterised on the 
fluid level in terms of a scalar resistivity. We have 
no reason to doubt this but of course a proof is be- 
yond the scope of these, or any currently available, 
numerical simulations. 

1. The instability will trigger the formation of an 
intense current concentration (a current sheet). 
The combination of instability and current sheet 
dissipation will cause the loop temperature to 
increase along its whole length by a factor of 3 
over the initial background value, i.e. heating 
up to about 6 x 10 6 K. A higher temperature 
component, perhaps 6 times the initial temper- 
ature (around 1.8 x 10 7 K), may also be observ- 
able (see Figure |A8|(b)). 

2. The whole resistive phase takes about 10 Alfvcn 
transit times (approximately 5 seconds) so the 
loop would brighten very rapidly. The viscous 
dissipation timescale for this loop is t„ ~ 4 min- 
utes. In this we have taken t v = L 2 v /v where L v 
is the velocity scale length and we have taken 
L v = 10 6 m. The conductive timescale (the 
timescale for which the loop should be visible 
at 3 times background temperature) is also of 
the order of minutes. 

3. For this equilibrium the total magnetic energy 
released is ~ 5 x 10 28 ergs. This is split approxi- 
mately equally between Ohmic heating, viscous 
heating and kinetic energy. The energy released 
for other size loops can be found from noting 
that the energy released scales as B^a\L z where 
Bq is the magnetic field in the loop, clq is the 
loop radius and L z the loop length. 



4. The loop is not destroyed but remains confined 
within the same region of the corona. The 
twisted field lines become straightened within 
the same confining region. 

5. After brightening there will be very large, lo- 
calised flows in the loop. These have a max- 
imum of ~ 0.8Va (~ 1500fcm s" 1 ). However, 
after taking into account the line of sight effect; 
exposure times and averaging over a diagnostic 
pixel size we find that the predicted observable 
flows are ~ 40km s -1 . The timescale of viscous 
dissipation means that these flows would per- 
sist for minutes after the initial brightening and 
therefore should be observable. 

6. The kinetic energy generated is about half the 
total thermal energy released (both Ohmic and 
viscous heating) and this may be indirectly ob- 
servable. The averaging inherent in measure- 
ments outlined above would mean that temper- 
ature estimates based on spectral line widths 
(which would included the averaged Dopplcr 
broadening) should be about 1.5 times those 
from calculations based on ratios of line inten- 
sities (which will be insensitive to the plasma 
motion predicted here). 

The clearest study of transient brightening loops 
(Shimizu ct. al. 1994) shows that the above diagnos- 
tic signatures are consistent with some of the SXT 
observations from Yohkoh. The energy released by 
the instability is at the high end of observed energies 
for brightening loops, i.e. they are micro-flares. This 
is as should be expected since the idealised nature 
of these simulations is only applicable to loops which 
have been twisted by a long timescale input of en- 
ergy. The majority of these loops are seen to brighten 
along thei r en tire length consistent with the picture 
in Figure ^8|(a). Assuming a background active re- 
gion coro nal temperature of 2 x 10 6 K the iso-surfaces 
in Figure A8(a) would correspond to temperatures of 
~6x 10 6 K . The more structured, high temperature 
component shown in Figure ^8|(b) would correspond 
to a temperature of ~ 2 x 10 7 K. The conductive loss 
timescale for a coronal loop at ~ 10 6 K is ~ 10 2 s, 
so one would expect the brightening of the loops pre- 
dicted here to last several minutes. However the very 
high temperature component is highly localised (see 
Figure A8) and would be smoothed out on a much 
more rapid timescale. It would therefore require an 
exposure time of the order of seconds, on a diagnostic 
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sensitive to temperatures ~ 2 x 10 7 K to confirm or 
dismiss the predictions shown in Figure |A8| (b) . The 
speeds predicted in these simulations are consistent 
with observed speeds from Doppler measurements but 
a more directed study is needed to confirm if the 



with 



size and structure shown in Figure A9 are present 
after a loop brightening event. The predicted dis- 
crepancy between temperature measurements based 
on line broadening and ratios of different spectral line 
intensities is also as yet untested. 

In summary, we have performed a set of numeri- 
cal simulations of unstable coronal loops which carry 
no net current. By using high resolution numerical 
experiments, with a resistivity given by equation |^, 
we have been able to have greater confidence in our 
quantitative predictions from these simulations than 
would have been possible by simply assuming uni- 
form resistivity. Our predictions for these observa- 
tional signatures are listed above. Where comparison 
with current observational data is possible these are 
in broad agreement with the properties of brightening 
(compact flare) loops. The full set of predictions can 
now be used as the basis of a directed observational 
study which may then confirm, or dismiss, whether 
large scale MHD instabilities are the cause of the high 
energy, micro-flare end of the spectrum of transient 
loop brightcnings. 
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A. Details of the Lagrangian Code 

The Lagrangian code used to follow the ideal evo- 
lution of the instability is based on the equilibrium 
code described in Longbottom et. al. 199£ and Craig 
It solves the ideal Lagrangian MHD 



k Sneyd 1986 



equations 



Dx 
~Dt 



(V x B) x B - VP, 



(Al) 
(A2) 



B, 
A 



constant (moving with the fluid)(,A3) 

(A4) 
(A5) 



Po/A, 
" Xi -B 0j /A 



dXj 
d(xi 1 x 2 ,x 3 ) 

d(x u x 2 ,x 3 y 



(A6) 



Here x = (xi, x 2 , x 3 ) is the current position of the 
fluid element which is initially at X = (Xi,X 2 , X 3 ), v 
is the velocity moving with the fluid, p and po are the 
current and initial densities of the fluid clement, B = 
(Bi 1 B 2 ,B 3 ) and B = (B i, B Q2 , B a3 ) arc the current 
and initial magnetic fields and P is the pressure. 



Equations Al and A2 are advanced in time using 



a Lax-Wendroff type method (fourth order in space, 
second order in time). Once the new positions of the 
fluid elements are known the other variables can be 



calculated directly from equations A3- A6 without fur- 
ther time integration. The method preserves total 
mass, entropy and V • B = identically and gives 
excellent energy conservation. 

As time progresses the grid deforms moving more 
points into regions of compression. For the time evo- 
lution considered here, with the inner part of the 
twisted loop being forced by the instability into the 
almost static external potential field, the grid points 
accumulate in the regions where the current sheet 
forms. This allows these increasing currents to be 
resolved for much larger values than would be pos- 
sible by the equivalent Eulerian code. It should be 
noted, however, that as the Lagrangian code relies 
on the system being ideal it can say nothing about 
the evolution of the system once dissipation becomes 
important. 
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Fig. Al. — Magnetic field components Bg and B z vs. 
radius for the initial equilibrium. 

Fig. A2. — The modulus of the current perturbation, 
|ji| at y = z = 0, t = 5, from the linear (dotted line), 
the non-linear Lagrangian (solid line) and non-linear 
Eulerian (dashed line) time-evolution simulations. 



Fig. A3.— As in Figure A2 but for t = 10 



Fig. A4. — The x-component of the velocity, v x at 
y = z = 0, t = 5, taken from the linear time evolution. 



Fig. A5. — The modulus of the current perturbation, 
|ji| at y — z — as a function of x and time from the 
nonlinear Lagrangian simulation. 

Fig. A6. — The modulus of the current, |j| at z = 
from the nonlinear Lagrangian simulation at t = 11.5. 



Fig. A7. — Iso-surfaces of \j\ = 3 at t = 10, and 
t = 15 from the Eulerian simulation. 

Fig. A8. — Iso-surfaces of 3 times background energy 
and 6 times background at t — 20 from the Eulerian 
simulation. 

Fig. A9. — Plasma velocity after averaging over a 10 
seconds exposure time and a 1.5 Mm x 1.5 Mm square 
vs. distance along the loop. Photospheric footpoints 
are at ±5Mm. 



This 2-column preprint was prepared with the AAS IATjrjX 
macros v4.0. 
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